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ABSTRACT 

We investigate the general interaction between an eccentric planet and a coplanar 
debris disc of the same mass, using analytical theory and n-body simulations. 
Such an interaction could result from a planet-planet scattering or merging event. 
We show that when the planet mass is comparable to that of the disc, the former 
is often circularised with little change to its semimajor axis. The secular effect 
of such a planet can cause debris to apsidally anti -align with the planet’s orbit 
(the opposite of what may be naively expected), leading to the counter-intuitive 
result that a low-mass planet may clear a larger region of debris than a higher- 
mass body would. The interaction generally results in a double-ringed debris 
disc, which is comparable to those observed in HD 107146 and HD 92945. As an 
example we apply our results to HD 107146, and show that the disc’s morphology 
and surface brightness profile can be well-reproduced if the disc is interacting 
with an eccentric planet of comparable mass 10 — 100 Earth masses). This 
hypothetical planet had a pre-interaction semimajor axis of 30 or 40 au (similar 
to its present-day value) and an eccentricity of 0.4 or 0.5 (which would since have 
reduced to ~ 0.1). Thus the planet (if it exists) presently resides near the inner 
edge of the disc, rather than between the two debris peaks as may otherwise be 
expected. Finally we show that disc self-gravity can be important in this mass 
regime and, whilst it would not affect these results significantly, it should be 
considered when probing the interaction between a debris disc and a planet. 

Key words: planets and satellites: dynamical evolution and stability - planet- 
disc interactions - circumstellar matter - stars: individual: HD 107146 


1 INTRODUCTION 

Given the calm order of the Solar System today, where 
most planets and minor bodies occupy near circular and 
coplanar orbits, one could be forgiven for forgetting that 
planetary systems can be violent places. Indeed, our own 
system probably had a tumultuous youth; planets may 


have scattered off each other and collided (Hartmann & 


Davis 


et al. 


19751), switched places with one-another (Tsiganis 


20051) and ploughed into regions of debris (Walsh 


e t al.| 2011). This turbulent picture is also inferred from 
extrasolar planets; many of these objects are eccentric 


(Schneider et al. 2011), a possible hallmark of previous 


scattering events (Juric & Tremaine 2008). Some, such 
as Fomalhaut b, may also pass through regions of de¬ 
bris ( Kalas et al.|2013 Pearce, Wyatt &; Kennedy|2015 ). 
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Furthermore, major orbital evolution is required to ex¬ 
plain some classes of extrasolar planets, such as the Hot 
Jupiters ( Rasio &; Ford|p~996 Wu &; Lithwick|[2011 ). So 
if planet-planet scattering, mergers and dynamical insta¬ 
bilities could be the norm then it is pertinent to ask how 
planets affected by these processes interact with other 
bodies in the system, and whether we can use this infor¬ 
mation to probe their past or ongoing dynamical evolu¬ 
tion. 

In this paper we examine the general evolution of a 
system hosting a debris disc interacting with an equal- 
mass, coplanar, eccentric planet, assuming the planet’s 
eccentricity was rapidly driven up by one of the above 
processes. Our Solar System hosts two debris discs, the 
Asteroid and Kuiper Belts, and many extrasolar discs 
have been detected as infrared excesses in the spectra 
of stars (e.g. Rhee et al. 2007| Eiroa et al. 2013). So 
given that debris discs are reasonably common, it is likely 
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Equal-mass planet disc interactions 2 


that dynamically evolving planets interact with these 
structures in at least some systems. Instruments such as 
Spitzer, the HST and ALMA have resolved some extra¬ 


solar debris discs (e.g. Backman et al. 2009 Schneider 


et al. 2009 Dent et al. 2014), and the eccentricity or 


dumpiness of these discs can be used to infer the pres¬ 
ence of planets which would be otherwise undetectable. 
Hence one aim of this work is to characterise this inter¬ 
action in general, to ascertain the signatures an eccentric 
planet leaves on a disc which can then be compared to 
observations. This part of the paper is similar in its goals 


and methodology to our previous investigation (Pearce & 
Wyatt||2014 ), in which we considered only planets much 


more massive than the disc; the main difference now that 
the planet and disc are of equal mass is the disc’s ability 
to significantly alter the planet’s orbit, which can have 
major effects on the system evolution. 

Our second aim is to apply the general results to the 
debris disc of HD 107146, which has been resolved in both 


Williams et al.| 2004; |Carpenter et al.12005 

fCorder et al. 

2009; Hughes et al. 2011; Ricci et al. 2015 

). This disc is 
) and axisym- 

seen nearly face on, and is broad (~ 100 au 


metric with a ~ 30 au inner hole. The curiosity of this 
system is that the 1.25mm debris surface density profile 
appears to either first decrease and then increase with 
radius, or generally increase but with a gap at around 


80 au (Ricci et al. 2015). This differs from the surface 


density profiles of protoplanetary discs, which decrease 
with radius (e.g Andrews fe Williams||200 7); it has been 
suggested that the unusual profile could be the result of 
planetary perturbations. Given the large disc mass (pos¬ 
sibly ~ 100M® in total; Ric ci et al.|2015 ), this system is 
a prime application of our results. We will show that the 
strange disc morphology can be explained as the after- 
math of the interaction between an eccentric planet and 
a coplanar debris disc of the same mass. 

The layout of this paper is as follows. We examine the 
general outcomes of this interaction in Section [2] using a 
combination of theory and n-body simulations. We then 
apply these results to HD 107146 in Section [3] where we 
run simulations over a broad region of parameter space to 
replicate the disc structure and surface brightness profile 
observed with ALMA. We discuss the implications of the 
results in Section^] and conclude in Section [5] 


2 GENERAL INTERACTION OUTCOMES 

In this section we describe the general outcomes of the in¬ 
teraction between an eccentric planet and a comparable- 
mass, coplanar debris disc. We first summarise the dy¬ 
namical processes involved in the interaction in sections 
|2.1| to |T3l and describe their effects on the orbits of both 
the debris particles and the planet. Wherever possible we 
give quantitative predictions from dynamical theory. In 
Section |2.4| we present an n-body simulation as an ex¬ 
ample of the general outcomes of this interaction, and 
explain these in the context of the theoretical results. 


2.1 Secular effects on debris 


Secular interactions are long-term angular momentum 
exchanges between bodies, which can cause a particle’s 
eccentricity and orbital plane (but not semimajor axis) 
to evolve. As secular timescales are much longer than or¬ 
bital timescales then, providing the objects are not in 
a mean-motion resonance, each interacting body may be 
thought of as an extended “wire” of material in the shape 
of the object’s orbit, with a density at each point inversely 
proportional to the body’s velocity there (this approach 
to secular perturbations is known as Gauss averaging). 
The influence of other masses causes this wire to change 
shape and orientation. There exists no general analytic 
evaluation of this interaction (although semi-analytic so¬ 
lutions can be constructed in some cases; see |Beust et al.| 
2014), so a common approach is to derive an analyti¬ 


cal form complete up to second order in eccentricities 
and inclinations, and disregard all higher terms. This is 
valid for small eccentricities and inclinations, but intro¬ 
duces errors if larger values are considered; however, we 


showed in Pearce & Wyatt (2014) that this still produces 


qualitatively correct results for very large eccentricities 
so long as the inclinations are small. We now summarise 
the behaviour of a particle undergoing a secular inter¬ 
action with a planet according to second-order secular 
theory, and apply the results to the interaction between 
an eccentric planet and a debris disc. 

Second-order secular theory predicts that a body’s 
eccentricity e and longitude of pericentre vj are coupled 


(for details see Murray & Dermott 1999). Specifically, 
for a test particle undergoing secular perturbations from 
one or more massive bodies, these quantities satisfy the 
equations 


e cos vj — e p cos uj p + ef cos tuf , (1) 

e sin vj = e p sin vo p + ef sin tUf. (2) 

Here e p and zu p denote the “free” parameters of the par¬ 
ticle; e p is constant, and zu p increases linearly with time 
t such that 


vj p = At + p, (3) 

where A and ft are constants. The quantities ef and vjf, 
are the “forced” values, which depend on the parame¬ 
ters of the massive bodies in the system and may evolve 
in time. Hence the test particle moves around a circle 
on the ecostu, esintu plane, of radius e p and at a rate 
A. Meanwhile the centre of this circle also moves as the 
perturbing bodies evolve over time. 

We now consider a system comprised of a star of 
mass M*, a test particle at semimajor axis a and a planet 
of mass Mpit at semimajor axis a p it (where a > a p it). We 
also assume that some mechanism causes the planet’s 
pericentre to precess (d7 p i t ^ 0). We set both the test 
particle’s eccentricity and the planet’s longitude of peri¬ 
centre to zero at time zero for simplicity, i.e. ft = 7r, 
vji — zjpit, e p = jefo | (where efo is ef evaluated at t — 0). 
The forcing eccentricity is now 
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Figure 1. Evolution of a test particle’s orbit under the influence of a single, precessing planet, according to second-order secular 
theory. Shown are the cases where its orbit precesses more quickly than that of the planet (A > ffj p i t ) and more slowly (A < ff7 p i t ), 
for a planet of a given eccentricity. Left plot: the coupled evolution of eccentricity e and longitude of pericentre vj, where the 
large black circle denotes the planet. A particle of a given semimajor axis will move around one of the two paths in the direction 
indicated, depending on the sign of A — d7 p i t . Right plots: the corresponding physical areas swept out by the particles, in a frame 
instantaneously aligned with the planet’s orbit. The asterisk and thick red line denote the star and the planet’s orbit respectively, 
and the dashed lines show the extreme orbits of the particle. They grey region is the area swept out by the particle, resulting from 
the superposition of all orbits between the two extremes. Both plots show systems with identical parameters, but with differing 
signs of A — rijpit- Hence a particle which would never cross the planet’s orbit if A > vj p \t may do so if A < zo p \t- Also note that 
the planet is precessing, so the structures on the right hand plots will rotate whilst remaining aligned with the planet’s orbit. 


e f = 




A — d 


^pit 


e pit? 


where 


ApH --- 


GM* Mpit Q-pit 7 ( 2 ) / / \ 

~ nr - m7^t 6 3 /2 («ph/«) 


(4) 


(5) 


1 / GM * M p it Q-pit (i) , / x 

A= 4\ < 6 > 

G is the gravitational constant and bi j \a) is a Laplace 
coefficient, such that 


1 r 2 

^ j \a ) = - / 
n Jo 


cos(j , 0)d'0 


(7) 

(1 —2a cos ^a 2 ) 8 v y 

Equation [4] shows that the orbit of the test parti¬ 
cle will evolve differently in the regimes A > ff7 p i t and 
A < zu pit. Firstly, if A > ffjpit then the precession rate of 
the particle is faster than that of the planet, and ef > 0. 
The particle will therefore move anticlockwise about a 
circle on the ecos^ — zj p h), esin(tu — tu p i t ) plane, which 
crosses the origin and is offset from the origin in the di¬ 
rection of the planet. This is shown by the solid line on 
the left hand plot of Figure [l] Hence the particle’s eccen¬ 
tricity will be maximised when its orbit is aligned with 
that of the planet, and small when antialigned. Super¬ 
imposing all intermediate orbits shows that the particle 
will sweep out a broad, eccentric disc aligned with the 
planet’s orbit (central plot of Figure [l]). This allows the 
particle to attain high eccentricities and yet be shielded 
from scattering by the planet. For a complete summary 
of particle evolution in this regime, see Pearc e~fe Wyatt | 
(2014) and Faramaz et al. (2014). 


If A < ffjpit, the planet precesses more quickly than 
the particle and the forcing eccentricity becomes nega¬ 
tive. The particle will now move clockwise about a cir¬ 
cle on the ecos(zj — vj v \t), esm(zj — vj p \t) plane, again 
crossing the origin but offset in the direction opposing the 
planet. This is the dashed line on the left hand plot of Fig¬ 
ure [l] The particle’s eccentricity will now be maximised 
when its orbit is antialigned with that of the planet, and 
superimposing all intermediate orbits results in a broad, 
eccentric disc antialigned with the planet’s orbit (right 
hand plot of Figure [I]). Hence a particle which would 
never cross the planet’s orbit if A > U7 p i t may now do so, 
and could therefore be ejected by the latter. 

Equation [6] shows that the particle’s precession rate 
A is set by the perturber’s mass and the semimajor axes 
of the particle and perturber, all of which are constants in 
the secular problem. On Figure [2] we plot A as a function 
of test particle semimajor axis for an example system, 
which contains a precessing planet. Generally, A 0 if 
the test particle’s semimajor axis is very small or very 
large, and A —> oo if the semimajor axes of the particle 
and planet are similar. Hence if zo v \t ^ 0 then there will 
always be a region in particle semimajor axis space where 
A > ffjpit, and two regions where A < In this case 

a sufficiently broad disc of test particles will cover both 
the A > d7 p it and A < ff7 p i t regimes. Figure [3] shows a 
schematic of the resulting debris structure, for particles 
external to the planet. The result is a superposition of 
the debris structures on Figure [l] with the innermost 
particles in the A > ff7 p i t regime and the outermost in the 
A < ffjpit regime. There exists a crescent-shaped region 
devoid of debris, and also a location where A > d7 p i t and 
A < ff7 p it particles overlap. However this overlap does 
not necessarily correspond to debris overdensity, and it 
is unlikely to be a site of increased dust production; the 
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Figure 2. Precession rate of a test particle as a function of 
its semimajor axis, derived using second-order secular theory. 
Here the planet (black circle) has a semimajor axis of 40 au 
and precesses at a rate dj p it = 4 x 10 -5 °yr _1 . Particles pre- 
cessing more rapidly than the planet form an eccentric disc 
apsidally aligned with the planet’s orbit (central plot on Fig¬ 
ure^), whilst those with slower precession rates will antialign 
with the planet’s orbit (right hand plot on Figure[l]). Note that 
this plot does not take account of mean motion resonances, the 
effect of which could dominate over secular behaviour. 



Figure 3. General shape of a debris disc undergoing a secular 
interaction with a interior, precessing planet. The grey regions 
are disc material, the red ellipse denotes the planet’s orbit and 
the asterisk shows the star. There are two distinct debris popu¬ 
lations; an inner population apsidally aligned with the planet, 
and an outer one which is antialigned. The inner population 
may or may not be present depending on the parameters of 
the system, and likewise for the outer population. 


collision velocities between particles are actually greater 
within the inner (aligned) disc than those in this inner- 
outer disc overlap region. 

One mechanism which may cause planet precession 
is the secular effect of the debris. If the planet is much 
more massive than the total disc mass, the precession 
rate of the former will be slower than that of the de¬ 
bris, unless the disc extends very close to or far from the 


star. Hence for all but the broadest discs, if the planet is 
much more massive than the disc then the result will be 
an eccentric debris structure apsidally aligned with the 
planet’s orbit. This was the case investigated in |Pearce| 


& Wyatt (2014). Alternatively if the planet mass is com¬ 


parable to that of the disc (ignoring planet evolution for 
now), the planet will likely precess more rapidly than 
the outermost debris. Thus the farthest debris will as¬ 
sume an eccentric structure antialigned with the planet’s 
pericentre. Whether the innermost debris also forms this 
structure, or forms a structure apsidally aligned with the 
planet (as on Figure [3]), depends on the parameters of 
the system; a sufficiently eccentric planet will eject all 
particles with similar semimajor axes, leaving only dis¬ 
tant debris which may be slowly precessing. In this case 
only the antialigned debris structure would be formed. 
This also leads to a counter-intuitive result; a planet of 
comparable mass to the disc will clear a larger region 
of debris than a much more massive planet. This is be¬ 
cause particle orbits antialign with that of a low mass 
(i.e. rapidly precessing) planet, and are therefore more 
likely to be ejected than if under the influence of a more 
massive planet (which their orbits would align with). 

Thus far we have only considered the secular evolu¬ 
tion of debris which does not intersect the planet’s or¬ 
bit. Particles with orbits crossing that of the planet will 
eventually be scattered unless in a mean-motion reso¬ 
nance, however secular evolution may still occur before 
then. Beust et al. (2014) simulated the evolution of de¬ 
bris under the influence of an eccentric planet, when the 
planet’s orbit crosses that of the debris. They showed 
that the secular interaction still drives up particle eccen¬ 
tricities as before. However their orbits do not preferen¬ 
tially align or antialign with that of the planet, but rather 
initially orientate themselves such that for much of the 
time their pericentres are misaligned with the planet’s 
by ~ 70°. We observed similar evolution of particles on 
planet-crossing orbits in our simulations using high mass 
planets (Pearce & Wyatt 2014), suggesting that this be¬ 
haviour arises because the orbits intersect. That we are 
now concerned with comparable planet and disc masses 
does not make a difference; the particles affected by this 
mechanism have semimajor axes similar to the planet’s, 
and therefore still precess faster than the latter even if 
the planet is rapidly precessing. Thus in addition to the 
long-term secular structures described above, particles 
crossing the planet’s orbit will have their eccentricities 
driven up and their orbits misaligned with the planet’s 
by ~ 70°, before eventually being scattered. 


2.2 The effect of scattering on debris 

Material which regularly crosses the planet’s orbit (again, 
if not in a mean-motion resonance) may be scattered by 
the planet. A particle’s post-scattering orbit may differ 
significantly from its pre-scattering orbit, but both must 
pass through the scattering point. Hence a planet which 
scatters material at all points around its orbit will form 
an overdensity of debris tracing its orbit, caused by the 
overlapping orbits of scattered particles. This overdensity 
is often strongest at planet apocentre, where scattering 
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is most efficient; this is because the planet spends more 
time around apocentre than at other points in its orbit, 
and the relative velocity between the planet and (non- 
scattered) disc particles is smallest here. 

Objects repeatedly scattered by the planet will even¬ 
tually leave the system or collide with other bodies. In 
the meantime, scattered objects may attain very large 
semimajor axes and eccentricities, and hence their orbits 
could extend far beyond the initial outer edge of the disc. 
Inclinations are typically excited less than eccentricities 
in repeated scattering encounters ( Ida &; Makino|p~992 ), 
so scattered material will form a broad disc superimposed 
on the secular debris structure discussed in Section [2711 

The surface density E of such a scattered disc will 
follow an r -3 ' 5 profile. This is an empirical result which is 
observed in all of our simulations regardless of parameters 
(as well as those of Duncan, Quinn &; Tremaine 1987), but 
it can also be obtained using the following semi-analytic 


method. According to the model of Yabushita (1980), 


particles repeatedly scattered by a planet will diffuse in 
semimajor axis space, such that the number of particles 
n with x in the range x x + dx (where x = 1/a) at 
time t is given by 


(l+ \/x/xo ) 


— \/x/xo 

T 


n(x, t) — — exp 

7 \ / 7 

l8) 

Here xo is the initial value of x, I 2 is the modified Bessel 
function, and r = t/tn{x 0 ) where to(x) is the diffusion 
timescale: 


t D (x) = 0.01T pltV ^x (^f) , (9) 

where T p i t is the orbital period of the perturbing planet. 
A reasonable approximation is that scattered particles 


tern, whereby debris particles initially on circular orbits 
with semimajor axes equal to a p it diffuse in x, whilst their 
pericentres remain at a p i t . We calculate the distribution 
of x values at time t, where t £r>(x 0 ), using Equation 
[8] We then create a virtual scattered disc consisting of 
a large number of particles, with semimajor axes drawn 
from the above distribution and pericentres equal to a p it. 
For each orbit we calculate the instantaneous radial dis¬ 
tance r of the particle at a randomised mean anomaly, 
and calculate the surface density profile resulting from 
the summation of these r values for all particles. Regard¬ 
less of the parameters used (planet mass, semimajor axis 
and stellar mass) this profile always goes as r -3 ' 5 . Hence 
an r -3 ' 5 profile appears to be a natural consequence of 
scattering, and a population of scattered material could 
potentially be identified from such a slope. 


diffuse m x whilst their pericentre distance and incli¬ 


nation remain constant (iDuncan, Quinn & Tremaine 


1987). Hence we may build a simple model of the sys- 


2.3 Planet evolution 

Unlike when the planet is much more massive than the 
disc, if the two are of comparable mass then the planet’s 
orbit may undergo significant evolution. It was noted in 


Section [2TT| that secular perturbations from the disc cause 
the planet’s orbit to precess, and its orbital plane will 
also evolve if the planet and disc planes are initially mis¬ 
aligned (although no significant plane evolution will occur 
if the two are roughly coplanar at the start of the inter¬ 
action). However the planet’s eccentricity may evolve sig¬ 
nificantly, through planet-particle scattering and secular 
interactions with the disc. These mechanisms individu¬ 
ally affect eccentricity in different ways, so overall the 
planet’s eccentricity behaviour combines two effects. Sec¬ 
ular perturbations from the disc will cause the planet’s 
eccentricity to increase and decrease periodically, whilst 
scattering damps the eccentricity and will circularise the 
orbit (given enough scattering events). Hence the planet’s 
eccentricity will undergo a long-term decline, with ad¬ 
ditional oscillatory behaviour in the meantime. If the 
planet scatters sufficient material before circularisation 
then there may be too little debris remaining to continue 
the damping process; in this case, the planet’s eccentric¬ 
ity may not tend to zero but to some higher value. 

Scattering will also change the planet’s semimajor 
axis. For a single planet scattering debris, the lack of inte¬ 
rior planets to remove material scattered inwards means 
than particles may only leave the system through colli¬ 
sions or ejection. The former mechanism will be rare as 
the star and planet pose small targets, hence the even¬ 
tual location of scattered material is likely exterior to 
its initial orbit. The planet will lose energy to counter 
this increase in particle energy, hence its semimajor axis 
will tend to decrease. In section 5.3 of Pear ce &; Wy-| 
att (2014) we derived a theoretical upper limit on this 


semimajor axis change, and showed that a planet cannot 
undergo significant migration if much more massive than 
the disc except for a contrived set of circumstances. The 
same arguments still apply even when the planet and 
disc are of comparable mass; in the context of the pa¬ 
rameters in equations 14-16 of Pearce &; Wyatt (2014), 
we require T ~ 1 for significant migration, which is un¬ 
likely for broad discs. Hence scattering is unlikely to cause 
any significant change in planet semimajor axis. Recall¬ 
ing that secular interactions also have no effect on this 
quantity, we conclude that the semimajor axis of an ec¬ 
centric planet interacting with a comparable-mass debris 
disc is unlikely to evolve significantly. 


2.4 General numerical simulations 


We now present n-body simulations of an eccentric planet 
interacting with a comparable-mass coplanar debris disc, 
to demonstrate the physical effects described in Sections 
2.1| to [273] We ran almost 100 n-body simulations using 
the Mercury 6.2 integrator (Chamb ers||1999 ), covering a 
broad region of parameter space. The general simulation 
setup is as follows. A planet of mass M p i t orbits a star of 
mass M*, with an initial semimajor axis a p it and eccen¬ 
tricity e p i t . The planet’s pericentre is typically of order 
1-10 au; this is roughly the location of the water snow 
line for solar-type stars, and hence the region where giant 
planets may be expected to form. The planets have initial 
eccentricities ranging from 0.1 to 0.9. We fix M* = 1M©; 
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changing this parameter affects the timescales in the in¬ 
teraction, but not the nature of the evolution. 

The star also hosts a debris disc exterior to the 
planet’s pericentre, of mass Mdi sc (where M p it = Mdi sc ), 
composed of N equal-mass particles. The disc midplane 
lies in the planet’s orbital plane. We consider discs with 
initial inner and outer radii (n and r<i respectively) of 
the order of 10 — 100 au. The semimajor axes a of disc 
particles have initial values between r± and 7 * 2 , and are 
distributed such that 


i(a) 


1 —7 

nr ' 


( 10 ) 


where 7 is the surface density index. These particles are 
initially on circular orbits, which are randomised in lon¬ 
gitude of ascending node and have inclinations up to an 
opening angle I with respect to the disc midplane. We 
use a pre-interaction opening angle of I = 5°, that of the 
classical Kuiper Belt ( Bernstein et al.|2004 ), and 7 = 1.5, 
that of the Minimum Mass Solar Nebula ( Hayashi|1981 ), 
in our simulations. 

We probe disc and planet masses from 0.1 Earth 
masses (0.1 M®) to 3 Jupiter masses. The discs contain 
N = 10 3 — 10 4 equal-mass debris particles, representing 
the more massive bodies (i.e. those unaffected by radi¬ 
ation pressure and PR-drag); hence only gravitational 
forces are included. Each particle exerts a force on the 
planet and vice-versa , but does not perturb other debris. 
Thus we ignore the self-gravity of the disc; this is dis¬ 
cussed in Section T4.1 1 Each simulation lasts of the order 
of 10 - 100 Myr. 

Despite the broad range of parameters tested, the 
interaction always produced the same qualitative results. 
We found four distinct evolutionary stages occurring on 
logarithmic timescales, which we describe below. We also 
present an example simulation; we show the disc surface 
density at each evolutionary stage on Figure [4] and the 
planet’s eccentricity evolution on Figure [5] The example 
shows a 10 M® planet interacting with a comparable- 
mass disc, with a p i t = 40 au, e p i t = 0.6, r± — 50 au, 
r 2 = 150 au and N = 10 3 . 


Stage 1: The planet begins to scatter material from 
the inner regions of the disc, depleting the debris surface 
density inwards of planet apocentre. Non-scattered par¬ 
ticles with orbits crossing that of the planet have their 
eccentricity increased via the planet’s secular influence; 
these orbits are preferentially misaligned by ± ~ 70° to 
the planet’s orbit, due to the effect described in |Beust| 
et al. (2014). Particles beyond the planet’s apocentre 


still have roughly circular orbits, and the similar secular 
phases of neighbouring particles cause the formation of 
a spiral-shaped overdensity beyond the planet’s orbit 


(Wyatt 2005). The planet’s eccentricity undergoes its 
most rapid decline due to debris scattering, and secular 
effects may also cause this eccentricity to oscillate. 


Stage 2 : All debris initially crossing the planet’s orbit 
has been scattered at least once; an overdensity of 
scattered material forms along the planet’s orbit, and 
this overdensity is strongest around planet apocentre. 
A population of scattered material with surface density 



z 


z 
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I 

CD 
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Figure 4. Example n-body simulation of an interaction be¬ 
tween an eccentric planet and an equal mass, coplanar debris 
disc, with the simulation parameters described in the text. The 
left panels show the debris surface density and the planet’s 
orbit (white ellipse), at time zero and then at subsequent evo¬ 
lutionary stages. The right panels show the radially averaged 
surface density at these times; the thick black lines are the sur¬ 
face density profiles, the dashed lines are the analytic surface 
density at t = 0, and the points show the planet’s semimajor 
axis and pericentre/apocentre distances. The red line on the 
Stage 3 and 4 surface density plots shows an r -3,5 profile, 
typical of scattered debris, and material beyond 150 au fol¬ 
lows this profile. The evolutionary stages are common to all 
our simulations, and are described in Section |2.4| 
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Figure 5. Eccentricity evolution of the planet in the simula¬ 
tion shown on Figure|4| as an example of the general behaviour 
observed in all our simulations. At early times (up to ~ 10 
Myr) secular eccentricity oscillations are noticeable, on top of 
the long-term decline from debris scattering. The dotted lines 
and numbers in boxes refer to the stages of system evolution, 
for comparison with Figure^] 


orbit, extending beyond the initial outer edge of the 
disc. This population is initially small, and hence is not 
visible on Figure [4] until the final panel. However a loga¬ 
rithmic surface density plot demonstrates that a ~ r -3 ' 5 
population has begun forming by the second stage. 
Material originating just exterior to the planet’s orbit 
may form a coherently eccentric disc apsidally aligned 
with the planet, depending on the planet’s precession 
rate and eccentricity (see Section 2 . 1 ). The spiral-shaped 
overdensity continues to develop beyond the planet’s 
orbit. The debris surface density profile hence has two 
peaks: a broad peak of scattered material stretching 
between the planet’s initial pericentre and apocentre 
distances, and sharp peak farther out corresponding to 
the spiral overdensity. 


Stage 3: At least one secular period has elapsed for 
particles initially orbiting just exterior to the planet’s 
apocentre; some initially stable material has been driven 
onto eccentric orbits apsidally antialigned with the 
planet, crossed its orbit and been scattered. Hence the 
surface density of scattered material tracing the planet’s 
orbit is increased, and a large crescent shaped gap 
forms in the disc in the direction of planet pericentre. 
The spiral overdensity exterior to the planet continues 
to move outwards, as more distant particles are still 
in secular phase with their neighbours. The planet’s 
eccentricity may by now have reduced significantly. 
Note that the planet’s location corresponds to a region 
of overdensity in the disc, rather than the region of 
underdensity (as might naively be expected). 

Stage 4: The planet has scattered all material crossing 
its orbit, and its eccentricity evolution essentially ceases. 


Hence the innermost peak of the surface density pro¬ 
file has been reduced or even removed. An overdensity 
of scattered material may still exist just exterior to the 
planet’s orbit; this material no longer comes close to the 
planet since the latter’s eccentricity decreased, so this de¬ 
bris is now stable. Particles driven to high eccentricities 
by secular effects early on may now have their eccen¬ 
tricities frozen, as the forcing eccentricity becomes small 
owing to the decrease in planet eccentricity. If the planet 
circularisation timescale is much longer than the secu¬ 
lar timescale of the outermost particles then surviving 
non-scattered debris forms a smooth disc apsidally an- 
tialigned with the planet; otherwise, the spiral overden¬ 
sity may still be present in the outer debris and remain 
there indefinitely. 


Generally, whilst the planet’s eccentricity and longi¬ 
tude of pericentre evolve significantly throughout the in¬ 
teraction, its other orbital elements remain roughly con¬ 
stant. In the example simulation a p i t changes by less than 
5 per cent, and z p it never exceeds 0 . 6 ° (from an initial 
value of 0 °). 

The qualitative results presented in this section are 
general. The quantitative results will differ for specific 
systems, but rough scaling rules can be applied. For ex¬ 
ample, increasing the mass of the planet and disc simul¬ 
taneously will decrease the interaction timescales, whilst 
increasing the planet semi-major axis and disc radii will 
increase timescales. Increasing the planet eccentricity will 
decrease the circularisation timescale, and moving the 
disc mass inwards (either through reducing rq or increas¬ 
ing 7 ) makes the planet circularise faster and to a greater 
degree. 

Whilst this paper primarily considers the case where 
Mpit = Mdisc, the results are applicable to the M p \ t > 
-Mdisc regime too. Even if the planet were orders of mag¬ 
nitude more massive than the disc, the general secular 
behaviour is the same as in the equal mass case. A differ¬ 
ence between the two mass regimes is that the transition 
between aligned and antialigned particles occurs farther 
from the star if the planet is more massive than the disc; 
this is because the planet would precess more slowly rel¬ 
ative to debris than the equal mass case, so the location 
where particles precess more slowly than the planet is far¬ 
ther from the star (see Figure [ 2 ]). This is why we did not 
observe this secular behaviour in Pearce & Wyatt (2014); 
our discs simply did not extend far enough outwards to 
probe this regime. The main qualitative evolutionary dif¬ 
ference between the equal mass case and that when the 
planet is much more massive is that the planet will not 
undergo the same degree of orbital evolution in the latter 
regime. 

An interesting result may occur if 1 < M p it /Mdisc < 
10 , whereby particles can change between the aligned 
and antialigned secular regimes. This occurs because the 
planet initially precesses rapidly (leading to antialign¬ 
ment of some particle orbits), yet the planet is massive 
enough to eject a significant fraction of the disc particles. 
The declining disc mass causes the planet precession to 
slow, meaning that the precession rate of some particles 
can “overtake” that of the planet. The final result of such 
an interaction is the formation of a coherently eccentric 
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disc (as in Pea rce fe Wyatt|2014 [, but with a more messy 
structure due to additional particles which have changed 
their secular behaviour. We do not wish to comment on 
the case where M p i t < Mdi sc , as disc self gravity would 
be very important in this regime and thus the results of 


this paper probably do not apply there (see Section 4.1) 


Our results may be used to predict the outcome of an 
eccentric planet interacting with a coplanar debris disc of 
the same or greater mass. They may also be used to infer 
the presence of an unseen perturber from the structure 
of an imaged debris disc, as we will now demonstrate for 
HD 107146. 


3 APPLICATION TO HD 107146 

HD 107146 is a 80-200 Myr old G2V star, located 27.5 pc 
from the Sun ( van Leeuwen|[2007 Williams et al.|[2 004). 
In 2000, IRAS imaging revealed excess infrared emission 


in the stellar spectrum, indicative of a debris disc (Silver- 
storie][2000). As noted in Section [l] further observations 


resolved the disc in both infrared emission and scattered 


light (Ardila et al. 2004 

Williams et al.|2004 Carpenter 

et al. 2005; Corder et al. 

2009 Hughes et al.|2011). For an 


excellent summary of work on HD 107146 up until 2011, 
see 


Ertel et al. (2011). 


The recent 1.25mm ALMA image reveals the disc 


at millimetre wavelengths in unprecedented detail (Ricci 
e t al.| 2015). These data show that the disc spans 30 — 150 


au from the star and, assuming it is circular, inclined by 
21° to the sky plane at a position angle (E of N) of 140°. 
These observations detected 0.2M® of dust at 1.25mm, 
and by extrapolating this up to bodies of diameter D — 
1000 km (with the number of bodies of a given diameter 
n(D) oc D~ 3 ' 6 ) the authors inferred a total disc mass 
of 100M®. The ALMA image and corresponding surface 
brightness profile are shown on the top two plots of Figure 

El 

It is clear from these plots that the disc has an un¬ 
usual morphology. The outer regions are brighter than 
the inner, and the brightness profile decreases with ra¬ 
dius before increasing again farther out. Lower-resolution 
880 /im SMA observations also show that the surface 
brightness does not decrease with radius as expected 
(Hugh es et al.p Oll). These data have been interpreted 
as the disc’s surface density profile either being double- 
peaked or increasing with radius with a gap at 80 au 


(Ricci et al. 2015), and these models are indistinguishable 
at the observation resolution. Possible causes of these pro¬ 
file include embedded Pluto-sized objects inducing colli¬ 
sions between large debris bodies, or perturbations from 
an unseen planetary companion. The ALMA observations 
failed to detect any CO gas, suggesting that a dust-gas 
interaction is not responsible for the disc morphology. 

We wish to ascertain whether a past (or ongoing) 
interaction between the disc and a hypothetical eccen¬ 
tric planet can explain the disc features, and if so, esti¬ 
mate the pre-interaction orbit of the planet as well as its 
present day location. We assume the planet originated 
interior to the disc, where some event placed it onto an 
eccentric orbit; this could have been a planet-planet scat¬ 
tering or merger event, for example (Lin & Ida 1997] Ford 


& Ra sio|2008 ). We aim to reproduce the 1.25mm ALMA 
observations with an n-body simulation of such an in¬ 
teraction. Millimetre grains are unaffected by radiation 
pressure and PR-drag, so should act as tracers of the par¬ 
ent debris bodies (those most important for the dynamics 
of the system). Hence the ALMA data is well-suited to 
modelling with purely gravitational n-body simulations. 


3.1 Simulation setup 

In addition to the simulations described in Section |2.4| 
we ran a further ~ 150 simulations specifically aimed 
at reproducing the HD 107146 debris disc. We describe 
their setup now. HD 107146 is a G2V star, so we fix its 
mass at 1M®. We also fix the disc mass to the 100M® 


value of Ricci et al. (2015). However this still leaves nine 


physical variables: the planet’s mass, its initial semimajor 
axis, eccentricity, inclination and argument of pericentre, 
and the disc’s initial inner and outer radii, opening an¬ 
gle and surface density profile. Computational limitations 
prevent us from exploring this whole parameter space, so 
we make several assumptions about the pre-interaction 
system to reduce the number of variables. 

We again fix the pre-interaction disc opening angle 
at I = 5° and 7 = 1.5, and again assume the planet 
initially orbits in the disc midplane. These assumptions 
leave five physical parameters: M p i t , a p it, e p it, r± and 7 * 2 . 
However we can use physical reasoning to fix a further 
two of these. Firstly, the disc of HD 107146 appears to 


be roughly axisymmetric. In Pearce & Wyatt (2014) we 
showed that if M p \ t ^ > -Mdi sc , the planet’s eccentricity 
will not be significantly damped by the disc, and external 
debris will form a coherently eccentric disc aligned with 
the planet’s orbit. Conversely, in Sections |2.3| and |2.4| we 
showed that if M p i t ~ Mdi sc then the planet’s eccentricity 
will be significantly damped, and hence the outer edge of 
the disc will remain roughly circular. Thus if HD 107146’s 
disc is interacting with a reasonably eccentric planet (or 
did so in the past) then M p i t ~ Mdi sc , so we fix the planet 
mass to be 100M® in our simulations. This means that 
the outer edge of the disc will be largely unchanged by the 
interaction, so we fix 7*2 = 151 au (which best fits the data 
in the outermost regions). Again, we use N = 10 3 — 10 4 
equal-mass debris particles to simulate the disc, and omit 
disc self-gravity (see Section 4.1). 

We are thus left with three free parameters: the ini¬ 
tial values of a p it and e p i t , and the initial inner disc radius 
7*1. These three are somewhat degenerate, so we cannot 
fix any at a single value. Instead, we run simulations with 
a p it = 20, 30, 40 and 50 au (noting that the inner peak of 
the observed surface density profile is at 50 au, and that 
the planet’s initial semimajor axis is typically interior to 
this peak in our general simulations), and for each a p it 
we run simulations with various values of e p it and n . We 
disfavour simulations where the planet’s initial pericentre 
is within 3 Hill radii of the disc inner edge, or external to 
this location; in these cases the disc would be unstable 
before the interaction started. Once our simulations are 
complete we compare them to the observations, using the 
method described below. 
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Figure 6. ALMA observations of HD 107146, along with our best-fitting simulation. Top left: ALMA 1.25mm continuum image 
( |Ricci et al. 2015J. Note that we use a different colour scale to that in the aforementioned paper. The white ellipse represents 
the beam size and orientation. Top right: points show the normalised, radially-averaged surface brightness profile of the disc as 
observed by ALMA, measured using elliptical apertures. The solid line is the profile from our best-fitting n-body simulation, at 
the time (19 Myr after the start of the interaction) of the best fit; the two agree with a reduced % 2 value of 0.4. The simulation 
parameters and the method used to compare the data and simulation are described in Section [3] Bottom left: positions of debris 
particles in the best-fitting n-body simulation at 19 Myr. The x — y plane is the initial disc midplane, with planet pericentre 
initially pointing along the x axis. The orbit of each particle has been populated with 100 points with randomised mean anomalies, 
to increase the effective number of particles plotted. The white point is the star, and the white ellipse the planet’s orbit. Bottom 
right: simulated ALMA image of the n-body disc. The particles have been scaled for emission, the image rotated, and smoothed 
with a 2D Gaussian representing the ALMA beam (white oval). Compare this to the ALMA observation in the top left, noting 
that we have not added noise and hence our image is smoother. 


3.2 Constructing simulated observations 


Throughout each simulation we compare the instanta¬ 
neous distribution of debris to that observed by |Ricci| 
(2015). This requires the simulated debris to be 


et al. 


converted into an image and surface brightness profile as 
would be observed by ALMA, for which we use the fol¬ 
lowing method. Firstly, we populate each particle’s orbit 
with 100 points at randomised mean anomalies, to in¬ 
crease the effective number of particles simulated. We 
then scale for emission by weighting each point by a 
black body; a point at radial distance r from the star 
is weighted to have a luminosity L, where 


L(r) oc B U (X, T). (11) 


Here B U (X,T) is the spectral radiance of a body of tem¬ 
perature T at a wavelength A, given by Planck’s law: 


B u ( A, T) oc 


exp 


he \ 
Xk B Tj 



( 12 ) 


where h is the Planck constant, c is the speed of light 
and k B is the Boltzmann constant. The temperature of 
the body is determined by the flux it receives from the 
star (again assuming black body behaviour), hence 


T =( ty ' 4 -~ i ' 2 < is) 

where L* is the star’s luminosity and a is the Stefan- 
Boltzmann constant. For HD 107146, we use L* — L® 
and A = 1.25mm, that of the ALMA observations. For 
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these parameters, B„{\,T) roughly scales as r~ 1 / 2 . For 
this analysis we have assumed the disc is optically thin; 
this is valid here since, whilst the disc is massive, it covers 
a broad region. Whilst the optical depth could be high if 
the disc were extremely thin, even a moderate 5° open¬ 
ing angle would be large enough that we do not have to 
consider flux attenuation here. 

To produce images for comparison with the ALMA 
observations, we rotate the simulated (emission scaled) 
disc to an inclination of 21° and a position angle of 
143°. We then convolve our image with a two-dimensional 
Gaussian to simulate the ALMA point spread function 
(PSF); this Gaussian has a standard deviation along its 
major axis of 13.4 au, along its minor axis of 9.8 au, and 
its major axis has a position angle of 19.8°. We also calcu¬ 
late the radially averaged surface brightness profile of the 
simulated disc, using elliptical apertures on the simulated 
image as in Ricci et al. (2015). We may then compare our 
simulations to the ALMA observations of HD 107146. 


3.3 Fitting the HD 107146 disc 


We identified the simulations which best replicate the HD 
107146 disc using a % 2 analysis. At 100 time intervals 
throughout each simulation we calculated the \ 2 value 
comparing the observed radial surface brightness profile 
with the simulated profile at this time (found using the 
method in Section 3.2). On Figure [ 7 ] we plot the a p i t , 
e p it, ri parameter space tested in our simulations, and 
colour each point by min(x? e d) (the minimum value of 
reduced x 2 > that is the minimum value of x 2 attained 
during that simulation, divided by the number of degrees 
of freedom). If the data are independent, a reduced x 2 of 
order 1 means that the obslervations are consistent with 
the model, and the smaller the value, the better the fit 
(although values much smaller than 1 imply the data is 
overfitted). Here the observed surface brightness profile 
points are correlated with each other, so little should be 
inferred from the exact value of min(xred) itself; however 
this value does allow a comparison between simulations, 
to identify that which best reproduces the HD 107146 
disc. 

Figure [ 7 ] shows that there are several regions of 
tested parameter space which produce discs consistent 
with ALMA observations. The best fit is attained using a 
planet with initial semimajor axis a p n = 40 au and eccen¬ 
tricity e p i t = 0.4, interacting with a disc with initial in¬ 
ner radius n = 50 au. For this case the simulated surface 
brightness profile is most similar to the ALMA observa¬ 
tions 19 Myr after the start of the interaction, and we plot 
the simulated ALMA image and surface brightness profile 
at this time on the lower two plots of Figure [6] These well 
reproduce the observations; the surface brightness profile 
yields a min(x? e d) value of 0.4, and the simulated image 
resembles the ALMA observation by eye. Note that glob¬ 
ular structures in the observed image are probably noise, 
which has not been accounted for in the simulated image 
and hence the latter appears smoother than the observa¬ 
tion. This best fit occurs when the simulated system is at 


stage 3 or 4 in its evolution (as described in Section 2.4); 
the planet has removed most of the material crossing its 


orbit, and its orbital evolution has essentially stalled. By 
this point the planet’s eccentricity has decreased from 0.4 
to 0.05, whilst its semimajor axis (initially 40 au) has only 
reduced by 4 au. The simulation first reaches x? e d ~ 1 at 
10 Myr, and this parameter remains less than 1 until the 
end of the simulation (at 30 Myr); hence the simulation 
also provides a good fit to the observations over a long 
time interval. 

However the best-fitting simulation is not unique in 
reproducing the observations. A well-defined x 2 mini¬ 
mum also exists for a planet semimajor axis of 30 au, cen¬ 
tred on e p it =0.55 and n = 50 au, and this is almost as 
good as the best-fitting 40 au solution (min(x? e d) — 0.5). 
Again, the planet in the a p i t = 30 au simulation un¬ 
dergoes minimal semimajor axis evolution whilst its ec¬ 
centricity is significantly reduced, and the simulation fits 
best once it has evolved to stage 3 or 4 and resembles 
the observations for a long time. Conversely, we find that 
planets with initial semimajor axes of 20 and 50 au do 
not reproduce the observed disc well. Also note that our 
well-fitting simulations have the disc’s initial inner edge 
exterior to its present day value, and material has since 
been scattered inwards by the planet. In conclusion, a 
planet with an initial semimajor axis of 30 or 40 au and 
an eccentricity of 0.4-0.5, interacting with a comparable- 
mass debris disc with initial inner edge at 50 au, can well 
reproduce the disc of HD 107146. At present, the planet 
is most likely on a roughly circular orbit at 30-40 au. 


4 DISCUSSION 

We have examined the general interaction between an 
eccentric planet and a coplanar, comparable-mass debris 
disc, and applied our results to HD 107146 in an attempt 
to explain its unusual disc. In this section we discuss a 
possible limitation of our work: the omission of disc self¬ 
gravity. We also discuss the timescale of the HD 107146 
interaction. Finally, we examine the implications of this 
paper for planet searches, both for general systems with 
debris discs and also for HD 107146. 


4.1 Disc self-gravity 


Our simulated debris particles exert a force on the planet 
(and vice-versa ), but do not interact with each other; 
hence we do not include disc self-gravity in our simu¬ 
lations. This omission dramatically increases computa¬ 
tional efficiency, allowing us to run several hundred sim¬ 
ulations for this paper. However whilst self-gravity does 
not affect the interaction outcome if the planet is much 
more massive than the disc (as in Pearce &; Wyatt|2014 ), 
if the two are of comparable mass then this effect could 
become important. 

Debris in a self-gravitating disc would undergo ad¬ 
ditional secular and scattering evolution from the influ¬ 
ence of other disc particles. Secular interactions work 
over large distances on timescales scaling inversely with 
the object masses, whilst scattering works over short dis¬ 
tances on timescales going as the inverse-square of the 
masses (see equations 17 and 18 in Pearce & Wyatt 2014). 
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Figure 7. Reduced x 2 of our simulated radially-averaged surface density profiles compared to that of HD 107146, at the time in 
each simulation when this parameter is minimised. We varied the initial inner disc radius r\ and the initial planet eccentricity e p i t 
for four initial planet semimajor axes a p i t , fixing all other parameters as described in the text. The points show our simulations, 
with the colourmap and contours interpolated between them. Contours show log 10 [min(x^ ed )] = 0, 0.5 and 1 respectively. The 
hatched regions show an unphysical area of parameter space, where the planet’s initial pericentre is closer than three Hill radii to 
the disc inner edge (see Section ED- 


Therefore given the small debris particle masses, the ma¬ 
jor effect of self-gravity is likely to be on the secular evo¬ 
lution of the disc. 

We investigate the possible secular effect of self¬ 
gravity by analytically calculating the precession rate of 
a test particle embedded in a disc. This gives us a feel 
for how the initial disc in our best-fitting HD 107146 
simulation would evolve due to self-gravity alone (in the 
absence of any planetary perturbations), and allows us to 
compare the magnitude of the self-gravity effect to that 
of the planet. 

To calculate the precession rate, we consider a test 
particle at position (R, 0, z) in cylindrical coordinates, 
which experiences a force from a 2 dimensional, axisym- 
metric disc in the z — 0 plane. Equation 2-146 in|Binney 


& Tremaine (1987) gives the radial acceleration of the 


particle due to the disc as 


ric disc exerts only a radial force on a disc particle. The 
precession rate of a particle at true anomaly / moving 
in a Keplerian potential and perturbed by an additional 
radial force F r is 


1 / a(l — e 2 ) 


F r COS / 


(16) 


(section 2.9 of Murray & Dermott |1999 ), and we average 
this over the orbital period T : 


W $ T l “ 2„,Wl-«» l ’ <17) 

where we have assumed that a and e are constant over 
one orbital period. Finally, if e 1 then F r will not vary 
significantly over the particle’s orbit. In this case 


F r (R) = ~ 


G 

R 3/2 


nr 

J r i 


K(k)~ 


1 k 2 


R R z' . , x 

X * ~R ~ ~Rf ~l~ JFJi ) ) 


4 1 -k 2 

kX(R')VR'dR f , (14) 


R R' ' R'R 
where R' is the radial location of a point in the disc, 


k 2 = 


4 RR' 


(R + R') 2 + z 2 ’ 


(15) 


and K(k) and E(k) are the complete elliptical integrals 
of the first and second kind respectively. We wish to con¬ 
sider a particle in the disc midplane, i.e. z = 0. However 
Binn e^T^ Tremaine| (1987) note that Equation |14| has an 
unphysical singularity at R — R' if z — 0, because k — 1 
here and so the K(k) and (1 — /c 2 ) -1 terms become unde¬ 
fined. This issue can be resolved by setting 0 < z « R, 
so we use z = 10 -4 au in our evaluation. A particle in 
the midplane experiences no vertical acceleration, and 
its tangential acceleration is also zero because the disc 
is axisymmetric. Hence the self-gravity of an axisymmet- 


(w> « (a), (18) 

V M 

and similar analyses for semimajor axis and eccentricity 
yield a « e ~ 0. We evaluated Equation |14| numerically 
for a disc with the initial parameters of that in our best¬ 
fitting HD 107146 simulation. The force, and the resulting 
precession rate, are shown as functions of radius by the 
black lines on Figure [8] 

The plot shows that the precession rate of particles 
inwards of 75 au is still dominated by the planet, even 
when disc self-gravity is considered. In the n-body sim¬ 
ulation (without self-gravity), the planet drives up the 
eccentricities of these particles and scatters the majority 
of them, with the remainder forming an eccentric disc 
aligned with the planet’s orbit. Hence this would still oc¬ 
cur with the inclusion of self-gravity. However the disc’s 
gravity may initially dominate beyond this region; in the 
n-body simulation, particles beyond 80 au preferentially 
antialign with the planet’s orbit, so debris out to 100 
au is removed. Self-gravity would effectively cause these 
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Figure 8. The evolution of a test particle under the influence 
of a massive axisymmetric disc. Top plot: the radial force im¬ 
parted by the disc on a coplanar particle at radius R , from 
Equation |14| The black line shows a 100 M® disc with an 
r ~ 1-5 surface density profile, with inner and outer radii of 50 
and 150 au respectively. These are the initial disc parameters 
in the best-fitting HD 107146 simulation. The red line shows 
the same disc but with all mass inwards of 75 au removed, 
representing the truncation of the disc by the planet. Bottom 
plot: the magnitudes of the resulting particle precession rates. 
The green dotted line shows the precession rate due to the sec¬ 
ular influence of a 100 M® planet with a p \ t = 40au (Equation 
[6j . The plot shows that the planet dominates particle evolu¬ 
tion in the inner regions of the disc, whilst disc self-gravity 
may be more important in the outer regions. 


particles’ precession rates to be uncorrelated with the 
planet’s evolution, preventing both preferential antialign¬ 
ment and also significant eccentricity excitation. So were 
disc self-gravity included, the depletion of the disc in the 
best-fitting HD 107146 simulation would initially extend 
out to 75 au rather than 100 au. This would not fit the 
observational data. 

However we have not yet considered the evolution 
of the disc self-gravity. Particles inwards of 75 au would 
be depleted even with self-gravity, and this would change 
the disc potential. The red lines on Figure [8] show a disc 
with the same parameters as discussed above, but with 
all mass inwards of 75 au removed. Now the planet is 
still influential out to about 95 au, so much of this debris 
may still eventually undergo scattering. Beyond this re¬ 
gion the disc self-gravity will always dominate, although 
in the simulation these particles were not significantly 
perturbed by the planet anyway. Hence the inclusion of 
self-gravity will not affect the overall simulation results 
in the outermost regions. 

Our analysis suggests that, for our best-fitting HD 


107146 simulation, the inclusion of disc self-gravity would 
not qualitatively affect the resultant disc structure inte¬ 
rior to 75 au and exterior to 95 au. In the region between 
these radii, the potential effect of self-gravity is unclear. 
This region might not undergo the same level of depletion 
as in the simulations, and the spiral structure visible in 
Figure [5] might not be present. Hence our observational 
fit might not be as good as that on Figure [6] More gener¬ 
ally, depending on the simulation parameters, self-gravity 
may affect our predicted outcomes for the interaction in¬ 
vestigated in this paper. The main effect of self-gravity 
would probably be the reduction of debris depletion in 
the region immediately interior to the outermost peak of 
the disc. However we stress that a more sophisticated self¬ 
gravity analysis is required to fully explore its potential 
effect, which is beyond the scope of this paper. 


4.2 HD 107146 interaction timescale 


Our best-fitting HD 107146 simulations all reproduce the 
observed disc ~ 10 Myr after the start of the interac¬ 
tion, compared to the 100 Myr age of the star. These 
two timescales are compatible, but scenarios in which 
they are comparable would be preferable. The interac¬ 
tion timescales are set by the disc (and hence planet) 
mass and, since the disc mass derived from observations 
is uncertain (Ricci et al. 2015), there is scope to change 
this in our simulations. The secular interactions between 
the planet and disc are the dominant effects in the sim¬ 
ulations, and Equations EMU show the secular preces¬ 
sion rate to scale linearly with mass whilst the forcing 
eccentricity is independent of mass. Equation [14] shows 
that the effect of disc self-gravity also scales linearly 
with disc mass, so scaling both M p i t and Mdi SC simul¬ 
taneously will not change the importance of self-gravity 
relative to planetary perturbations. Hence changing the 
disc and planet masses should affect the secular interac¬ 
tion timescales, but not the nature of this interaction. 
Reducing the masses will make the planet less efficient at 
ejecting debris, but seeing as the main effects of the in¬ 
teraction are secular in nature, this should not affect the 
outcome too much. Hence if we reduce the disc and planet 
masses in our simulations by an order of magnitude (so 
Mdisc ~ 10M®), then roughly the same interaction will 
occur over a timescale comparable with the stellar life¬ 
time. Hence whilst our interaction timescales are by no 
means incompatible with the system age, if we assume 
that this interaction is responsible for the observed disc 
structure then our results might suggest the disc mass 
is closer to 10M® than 100M®. Alternatively the planet 
may have only recently been placed on an eccentric orbit, 
and we happen to have observed the system at this stage 
in its evolution. 


4.3 Implications for planet searches 

Our findings have interesting implications for the infer¬ 
ence of unseen planets from debris disc features, both 
generally and for HD 107146. An important result is that 
the planets in this interaction generally circularise with 
little change in semimajor axis, having ejected much of 
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the debris interior to their final orbital distance. Hence 
the eventual location of the planet is typically near the 
inner edge of the disc, on an orbit traced by a debris over¬ 
density, beyond which lies a gap followed by another over- 
dense ring. This configuration would not otherwise be ex¬ 
pected; having observed a double-peaked debris disc, the 
nai've assumption would be that any perturbing planet 
lies in the underdense region between the two peaks. 
Hence future planet finding missions should not be dis¬ 
couraged if no planets are found in a debris disc gap; 
indeed, the absence of planets in this region may hint 
at a violent dynamical history, and could motivate the 
search for planets near the inner edge of the disc instead. 

If our hypothesis on the evolution of the HD 107146 
system is correct, then a ~ 100planet currently orbits 
near the inner edge of the disc, with a semimajor axis of 
30 — 40 au and an eccentricity of ~ 0.1. This planet orig¬ 
inated interior to the disc; assuming it was scattered out 
of its original location by another body then, based on its 
initial pericentre, a second companion with mass at least 
equal to that of the scattered planet exists at > 10 — 25 au 
from the star. Companions of less than 10 Jupiter masses 
(3OOOM0) have not been ruled out anywhere in the sys¬ 
tem by imaging (Apai et al. 2008), so this scenario is 
possible and could be tested with deeper planet searches. 
Furthermore, our simulations suggest that the outermost 
debris peak actually forms a thin spiral, rather than a 
continuous ring. This structure would be detectable in 
observations with ~ 3 times the resolution of the ALMA 
image, and such a detection would favour our hypothesis 
on the history of the system (although with the caveat 
that the disc self-gravity would also have an effect, and 
may partially or completely wash-out this spiral). Such a 
resolution may well be possible with current instrumen¬ 
tation. 

Another potential application of this work is to HD 
92945; this system may also harbour a double-peaked de¬ 
bris disc (Golimowski et al. 2011), so our results could 
be used to invoke a perturbing planet in that system 
too. However the HD 92945 disc was imaged in scattered 
light, so the emitting dust would be affected by radiation 
forces. Hence it is unclear without more detailed analysis 
whether the more massive debris also follows this double- 
peaked profile (as in HD 107146), or whether the observed 
morphology is a consequence of non-gravitational forces 
on small dust. 


5 CONCLUSIONS 

Broad, double-ringed debris discs could potentially have 
evolved to their present state under the influence of an 
eccentric, comparable-mass planet. We investigate this 
interaction in general, and show that it follows four dis¬ 
tinct stages on logarithmic timescales. A key result is 
that planet precession may cause distant debris orbits 
to anti -align with that of the planet, whilst the inner¬ 
most debris orbits align with the planet’s. This results in 
distinct inner and outer debris regions with a gap or de¬ 
pletion between them, akin to the double-peaked debris 
structures potentially observed in HD 107146 and HD 
92945. It also produces the counter-intuitive result that a 


low-mass planet may clear a larger region of debris than 
a higher-mass body. In general the planet undergoes a 
rapid eccentricity decrease whilst its semimajor axis re¬ 
mains constant; thus if the planet initially scattered off 
another body then the two would quickly decouple, so 
our results still hold in the presence of additional mas¬ 
sive planets (providing the eccentricity damping is fast 
enough). 

We then modelled the HD 107146 system in detail, 
confirming that the debris disc’s unusual morphology can 
be well explained by this interaction. If an unseen eccen¬ 
tric planet did sculpt debris into the structure seen today, 
then this hypothetical planet initially had pericentre in 
the inner regions of the system and apocentre within the 
disc itself; based on our best-fitting model, the planet 
is currently on a low-eccentricity orbit 30-40 au from the 
star. This is below the companion detection thresholds of 
current observations of the system, but could potentially 
be found by future imaging projects. 
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